/*=================================================================================
 * The functions of Mesh
 *=================================================================================*/

/*---------------------------------------------------------------------------------
 * Standard C++ library headers
 *---------------------------------------------------------------------------------*/

#include <fstream>
#include <iostream>

/*---------------------------------------------------------------------------------
 * Local headers
 *---------------------------------------------------------------------------------*/

#include "Mesh.h"
#include "FileReader.h"

using namespace std;

/*---------------------------------------------------------------------------------
	Name:		Mesh::MeshProcess

	Input(1):	The grid name

	Function:	Read the grid file to get the mesh;
				Initialization;
				Set the boundary;
				Check the mesh

	Return:		None
 *---------------------------------------------------------------------------------*/

void Mesh::MeshProcess(char *gridname)
{
	/* Read the grid file to get the mesh */
	FileReader grid;
	grid.readFile(gridname);
	ni = grid.getIntParameter("ni");
	xa = grid.getdoubleParameter("xa");
	xb = grid.getdoubleParameter("xb");
	nj = grid.getIntParameter("nj");
	ya = grid.getdoubleParameter("ya");
	yb = grid.getdoubleParameter("yb");
	bc = grid.getdoubleParameter("bc");

	/* Mesh initialization */
	Initial();

	/* Set the boundary */
	Boundary();

	/* Check the mesh */
	MeshCheck();
}

void Mesh::MeshProcess_MPI(int MPI_ID, int MPI_SIZE, char *gridname)
{
	int for_y_p = 0;
	/* Read the grid file to get the mesh */
	FileReader grid;
	grid.readFile(gridname);
	total_ni = grid.getIntParameter("ni");
	total_xa = grid.getdoubleParameter("xa");
	total_xb = grid.getdoubleParameter("xb");
	total_nj = grid.getIntParameter("nj");
	total_ya = grid.getdoubleParameter("ya");
	total_yb = grid.getdoubleParameter("yb");
	bc = grid.getdoubleParameter("bc");
	/* x-direction initiation*/
	ni = total_ni ;
	xa = total_xa ;
	xb = total_xb ;
	dx = (xb-xa)/(ni-1);
	xnode.Initial(ni+2*bc);
	total_xnode.Initial(total_ni+2*bc);
	for (int i = bc; i < ni + bc; i++)
		xnode(i) = xa + (i - bc) * dx;
	for (int i = bc; i < total_ni + bc; i++)
		total_xnode(i) = total_xa + (i - bc) * dx;

		
	/* y-direction initiation*/
	total_ynode.Initial(total_nj+2*bc);
	dy = (total_yb-total_ya)/(total_nj-1);
	for(int i=bc;i<total_nj+bc;i++)
		total_ynode(i)=total_ya + (i-bc)*dy ;
	int num_per_proc = floor(total_nj/MPI_SIZE);
	for_y_p = num_per_proc*MPI_ID ;
	ya = total_ya + for_y_p*dy ;
	if (MPI_ID == MPI_SIZE-1) num_per_proc = total_nj - (MPI_SIZE-1)*num_per_proc;
	nj = num_per_proc ;
	ynode.Initial(nj+2*bc);
	yb = ya +(nj-1)*dy ;
	for(int i=bc;i<nj+bc;i++)
	{
		ynode(i) = ya +(i-bc)*dy ;
		//total_ynode(for_y_p+i-bc) = ynode(i);
	}


	/* Set the boundary */
	Boundary();
	MeshCheck_MPI(MPI_ID);

	/* Check the mesh */
	//MPI_Gather(&total_ynode, ARRAY_SIZE, MPI_DOUBLE, total_ynode[0], , MPI_DOUBLE, 0, MPI_COMM_WORLD);
	if(MPI_ID ==0) TotMeshCheck();
	
	
	
}

/*---------------------------------------------------------------------------------
	Name:		Mesh::Initial

	Input(0):	None

	Function:	Initialize the mesh

	Return:		None
 *---------------------------------------------------------------------------------*/

void Mesh::Initial()
{
	xnode.Initial(ni + 2 * bc);
	dx = (xb - xa) / (ni - 1);
	for (int i = bc; i < ni + bc; i++)
		xnode(i) = xa + (i - bc) * dx;

	ynode.Initial(nj + 2 * bc);
	dy = (yb - ya) / (nj - 1);
	for (int j = bc; j < nj + bc; j++)
		ynode(j) = ya + (j - bc) * dy;
}

/*---------------------------------------------------------------------------------
	Name:		Mesh::Boundary

	Input(0):	None

	Function:	Set the mesh boundary

	Return:		None
 *---------------------------------------------------------------------------------*/

void Mesh::Boundary()
{
	for (int i = 0; i < bc; i++)
	{
		xnode(i) = 2 * xnode(bc) - xnode(2 * bc - i);
		xnode(ni + bc + i) = 2 * xnode(ni + bc - 1) - xnode(ni + bc - 1 - i - 1);
		ynode(i) = 2 * ynode(bc) - ynode(2 * bc - i);
		ynode(nj + bc + i) = 2 * ynode(nj + bc - 1) - ynode(nj + bc - 1 - i - 1);
	}
}

/*---------------------------------------------------------------------------------
	Name:		Mesh::MeshCheck

	Input(0):	None

	Function:	Check the mesh

	Return:		None
 *---------------------------------------------------------------------------------*/

void Mesh::MeshCheck()
{
	std::ofstream outfile("./output/grid_check.dat");
	outfile << "VARIABLEs = X,Y" << '\n';
	outfile << "ZONE I=" << ni + 2 * bc << '\t' << "J=" << nj + 2 * bc << '\n';
	outfile << "datapacking=block" << '\n';

	for (int j = 0; j < nj + 2 * bc; j++)
		for (int i = 0; i < ni + 2 * bc; i++)
			outfile << xnode(i) << '\n';

	for (int j = 0; j < nj + 2 * bc; j++)
		for (int i = 0; i < ni + 2 * bc; i++)
			outfile << ynode(j) << '\n';
}

void Mesh::MeshCheck_MPI(int id)
{
	std::ofstream outfile("./output/grid_check"+to_string(id)+".dat");
	outfile << "VARIABLEs = X,Y" << '\n';
	outfile << "ZONE I=" << ni + 2 * bc << '\t' << "J=" << nj + 2 * bc << '\n';
	outfile << "datapacking=block" << '\n';

	for (int j = 0; j < nj + 2 * bc; j++)
		for (int i = 0; i < ni + 2 * bc; i++)
			outfile << xnode(i) << '\n';

	for (int j = 0; j < nj + 2 * bc; j++)
		for (int i = 0; i < ni + 2 * bc; i++)
			outfile << ynode(j) << '\n';
}

void Mesh::TotMeshCheck()
{
	std::ofstream outfile("./output/grid_check.dat");
	outfile << "VARIABLEs = X,Y" << '\n';
	outfile << "ZONE I=" << total_ni + 2 * bc << '\t' << "J=" << total_nj + 2 * bc << '\n';
	outfile << "datapacking=block" << '\n';

	for (int j = 0; j < total_nj + 2 * bc; j++)
		for (int i = 0; i < total_ni + 2 * bc; i++)
			outfile << total_xnode(i) << '\n';

	for (int j = 0; j < total_nj + 2 * bc; j++)
		for (int i = 0; i < total_ni + 2 * bc; i++)
			outfile << total_ynode(j) << '\n';
}

/*---------------------------------------------------------------------------------
	Name:		Mesh::GetXnode/Ynode

	Input(0):	None

	Function:	Get the mesh nodes

	Return:		The mesh nodes
 *---------------------------------------------------------------------------------*/

Array<double> Mesh::GetXnode()
{
	return xnode;
}

Array<double> Mesh::GetYnode()
{
	return ynode;
}